Thermodynamics and dynamics of the formation of spherical lipidic vesicles 
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We propose a free energy expression accounting for the formation of spherical vesicles from 
planar lipidic membranes and derive a Fokker-Planck equation for the probability distribu- 
tion describing the dynamics of vesicle formation. We found that formation may occur as 
an activated process for small membranes and as a transport process for sufficiently large 
membranes. We give explicit expressions for the transition rates and the characteristic time 
of vesicle formation in terms of the relevant physical parameters. 

I. INTRODUCTION 

Phospholipid vesicles have been widely used as model systems for studying the dynamics and 
structural features of many cellular processes, such as endocytosis exocytosis Q], cell fusion 



[3, 4], transport and diffusion phenomena and membrane elastic properties j^]. 



In addition 



to its importance for basic research in the biological sciences, closed vesicles (liposomes) has been 
used as vehicles for the encapsulation of macromolecules such as nucleic acids {7, 



polymers and smal 
glass micropipette 



as well as 

molecules Large enough vesicles can be individually manipulated with a 



lOl . 111! ], and the vesicle membrane rigidity and, in general, membrane elastic 

properties can be measured 112 . Il3l . They have also been used as microreactors useful in the study 

I I I I 

of chemical reactions in geometrically confined spaces (141.1151]. In general, lipidic vesicles constitute 
nanocontainer systems ideally suited for the isolation, preservation, control and transport of a small 
number of molecules. 

There is a variety of experimental methods to prepare phospholipid vesicle suspensions, reviewed 
in [19J ■ One of the most widely used methods is the hydration of a dry phospholipid film 
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resulting spontaneously in a population of multilamellar vesicles with a high polydispersity in 
sizes and shapes. On the other hand, the formation of a unilamellar vesicle usually involves an 
intermediate structure in the form of a planar bilayer fragment, which is unstable, due to its exposed 



edges. These small planar bilayers can be grown 
or they can be formed from pre-existing bilayers 



jy detergent depletion, phospholipid precipitation 



la]. It is possible to prepare a population of giant 



unilamellar spherical vesicles when a dry phospholipid film is hydrated in the presence of an AC 
electric field jlij. I20I] . The resulting vesicle radius can be as high as 50 micrometers. A similar effect 
is exhibited by charged phospholipids. The bilayers ionize upon contact with water and they swell 
due to the repulsion of the bilayers, leading to the spontaneous formation of unilamellar vesicles 



2l(. 



In spite of the large experimental work existing, to our knowledge there is no systematic theo- 
retical model describing the dynamics of formation of a unilamellar spherical vesicle from a small 
planar membrane. Such a model could be useful for the characterization and control of the vesicle 
formation process and it could be tested by performing single- vesicle simulations and experiments. 
For instance, the video microscopy analysis of the closing dynamics of laser-generated transient 
pores on phospholipid membranes 



221 ] could be very useful in this regard. 
In this article, we propose a simple theoretical model for spherical-vesicle formation from a 
planar membrane, assuming that membrane rigidity and edge tension are the main contributions. 
We first calculate the free energy cost of vesicle formation and then, using this free energy and the 
rules of mesoscopic nonequilibrium thermodynamics (MNET), we derive a Fokker-Planck equation 
governing the evolution in time of a nonequilibrium distribution function that depends on time and 
the mesoscopic variable characterizing the instantaneous state of the system. Our analysis leads 
to identify that the ratio between the contour energy to curvature energy determines two main 
mechanisms of vesicle formation: i) An activated process for small values of the energy ratio and ii) 
a transport process for values larger than a critical value of the energy ratio. A detailed analysis of 
these two cases is performed leading to explicit relations for the vesicle formation rates in the first 
case and for the characteristic formation time in the second one. Our analysis is complemented 
with a numerical solution of the Fokker-Planck equation. 

MNET has been also used in other nanometric processes where curvature and surface tension 



33]. 



effects are the main driving forces, such as matter agglomeration systems; see for example 
The effect of linear tension on growth morphologies in 2D has been also studied in [jyj], where the 
entropy production has been shown to be the dominant selection mechanism. 

The article is organized as follows. Sec. [II] is devoted to derive the expression for the free energy 
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cost of vesicle formation by using equilibrium arguments. In section |III| we use this free energy 
to formulate a kinetic model for vesicle formation in terms of a Fokker-Planck equation and to 
analyze its implications. Finally, in section [IV] we present conclusions. 

II. THE FREE ENERGY 

In this section we formulate a simple model for the free energy associated to a phospholipid 
membrane in the process of wrapping in order to form a spherical vesicle. 

We will assume that in every stage of the process the membrane adopts the form of a spherical 
bowl as shown in Figure HJ In this process, we will consider two competing energies, one Fb 
associated with the bending of the membrane that favours planar membranes and another one Fi 
due to the contour of the membrane which favours spherical vesicles. 

According to the well established Helfrich theory the free energy of bending per unit area, 
obevi™ the nation F B = / fsdA, associated to a ioca! detonation ot a m e mbr ane is g ,ven by 

mm 

f B = 2K(H-c ) 2 + KK, (1) 

where k and K are the bending and the saddle-splay moduli respectively, cq is the spontaneous 
curvature of the bilayer. Here H = (l/2)(ci + C2) is the mean curvature, K = C1C2 is the Gaussian 
curvature and c\ and C2 are the local principal curvatures of the system. Since we are interested 
in homogeneous bilayers, then we may assume cq = 0. In our bowl approximation, both principal 
curvatures are identical and equal to inverse radius of the sphere r: c\ = C2 = 1/r. Therefore, the 
bending free energy simplifies to 

Fb = A^, (2) 

where Kb = 2k + k and A is the area of the membrane which will be assumed as constant. The 
contour free energy has the simple form 

Fi = Jl, (3) 

where 7 is the edge tension and I is the contour length. The total free energy is the sum of both 
contributions F = Fb + Eqs. ((3J and ([3|) can be rewritten in terms of the angle 6 (see Figure 
[TJ, leading to the following expression for the free energy 

F{6) = 2nK b [1 + cos(0)] + 2 7 (vr A) 1/2 sin(0/2). (4) 
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Planar membrane 



FIG. 1: Schematic representation of vesicle formation from a planar membrane. 



To derive this equation, we have used the fact that the total area of the bowl is A = 
2irr 2 [1 + cos(0)], that the contour length is given by I = 2nr sin(0), and used the trigonometric rela- 
tion sin(0)/ [1 + cos(0)] 1/2 = v / 2sin(6>/2). For convenience, we will use the following dimensionless 
form of the free energy 

F 



F(9) 



l + cos(0) + 5sin(0/2), 



(5) 



where S = {A/tt) 1 / 2 ^ / n b ). 

In Eq. ((SJ) it is clear that the parameter 5 determines the form of the free energy as a function 
of 0, and thus it determines when the planar membrane is stable and when it will spontaneously 
form a closed vesicle. We can identify the following three regimes (see Figure [2|): 

1. For 5 < 2, that is, when the linear tension is small compared to the bending constant 
then the free energy has a minimum at 9 = ir that corresponds to a planar membrane. 

2. For 2 < 5 < 4, there is a competition between contour and bending forces. As a result of 
this, the free energy has a minimum at 9 = corresponding to a closed spherical vesicle with 
an energy barrier centered at 9* = 2 arcsin(5/4). The free energy difference with respect to 
the planar membrane is given by 

1 



AF = F{9*) - F(ir) = -(5 - 4) z . 



(6) 



3. For 5 > 4 there is no energy barrier and the closed spherical vesicles are formed sponta- 
neously. 

Let us now estimate the possible values of 5 for real systems. For lipid bilayers the typical 
experimental values of the bending modulus are n ~ 5 — 25/cgT whereas for block copolymer 



bilayers a typical value is k ~ 40/csT, 
given by k ~ — an with a ~ 1 or less 



the order of 7 ~ 1 — 2 hsT/rim, 



25 



27, 



26] 



26j. The bilayer saddle-splay modulus is approximately 
2j|. Therefore, Kb ~ 5 — 25/cbT. The edge tension is of 



FIG. 2: Dimensionless free energy F(8) as function of 9 for 6 — 0.1, 1.8, 2.8, 4.1. 



For definiteness let us consider ~ 25ksT and 7 ~ lksT /nm leading to minimum radius of 
the vesicles (corresponding to 5 = 2) of r m j„, ~ = 25nm. For radius in the range between 
r ~ 25 — 50nm an energy barrier has to be overcome in order to form vesicles while for radius 
larger than 50nm the vesicles will form spontaneously. 



III. DYNAMICS OF VESICLE FORMATION 



At isothermal conditions, the free energy given in Eq. ([5]) can be interpreted (up to a constant) 
as the energetic cost or the minimum work necessary to form a vesicle 



Rr, 



2-TTKb [1 + cos(0) + <5sin(0/2)] . 



(7) 



This quantity can be used to derive a Fokker-Planck equation for the distribution function 
P(0, t) of finding the membrane in a stage characterized by 6 at time t. This distribution function 
is normalized and then satisfies a continuity equation of the form 

dP d 



(Pve), 



dt dd 

where Pvq is a diffusion probability current in #-space. 



(8) 



The Fo: 
namics 1 29 



kker- 
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lanck equation can be obtained by using the rules of nonequilibrium thermody- 



311 ] . This objective can be achieved by first calculating the entropy production of 
the system using the Gibbs entropy postulate [30} 



As = -k B I Pln-^-d6 



P 



0) 



eq 



In this equation As is the entropy change in the process of formation of the vesicle and the 
integration is carried out over all the range of values of 9 (from zero to ir). Here, P eq is the 



equilibrium reference distribution given by 



32] 



p p r „-Rmin/k B T_ p , -2wK b [l+cos(e)+Ssm(6/2)]/k B T 
1 eq — - r c — r e j 



(10) 
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where Po is a normalization factor. When 5 < 2, Eq. (jlOp is a very narrow function centered 
around 9 = tt, thus implying that the system remains as a planar membrane. For values of 5 
slightly larger than 2, one finds that the equilibrium state has a coexistence of vesicles and planar 
membranes. Otherwise, closed vesicles are the preferred configuration of the system. 

Now, by taking the derivative of Eq. ([9j) with respect to time and using Eq. ([8]), we obtain for 
the time derivative of the entropy: 

r^(pvein^-)de-i rpvMdo. (ii) 



dt D J de y p eq > tJ °de 

where we have defined the nonequilibrium chemical potential \i = fceTln \P/P eq \. This equation 
contains two terms, the first one constitutes the entropy flow and the second one is the entropy 
production a, given by a = (1/T) J Pvg(d/j,/d9)d8. From Eq. (jlip we may formulate linear 
relationships between the current Pvg and its conjugated force |^ in the form 

Pv 6 = -<*P%, (12) 



2SJ. This 



where a is the corresponding Onsager coefficient satisfying Onsager reciprocity relations 
use of a linear relationship assumes that the process is not too far from equilibrium, and therefore 
it may be not valid in general (this may be the case of charged vesicles). Note that for values 5 > 2, 
in the boundary 9 = 0, the velocity vg should vanish since the free energy has a local minimum. 
This condition does not affect the election (|12|) for other values of 9. 

Now, by substituting Eq. (fT2|) into ([8|), and using (flO|) we finally obtain 



dP d . 



sm(6>) + - cos I - 



89 

This is the Fokker-Planck equation governing the time evolution of the probability distribution 
during the formation of the vesicle. It contains a driving term characterized by the force 8F/d8 
and a diffusion term characterized by the diffusion coefficient D = ksTa. Here a plays the 
role similar to that of a friction or mobility coefficient in usual Brownian motion. In this case, 
it can be interpreted as a parameter characterizing the viscous or friction forces exerted on the 
membrane by the solvent or even it may include interlayer friction 35|. To make sure that P(9,t) 



will remain confined in the range < 9 < tt the initial condition can be written in the form 
P(9, 0) = [H{9) - H(tt - 9)]P (9) with H{9) the Heaviside function. 



A. Vesicle formation in the presence of energy barriers 

The Fokker-Planck equation (|13p can be used to calculate the transition rates r p from planar 
membranes to spherical vesicles (and for the inverse process r v ) in the regime where the energy 
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barrier controls the dynamics, that is, for values of 5 in the interval 5 £ (0,4). This objective can 
be achieved by following the usual methods of activated processes [36j. 

Let us estimate the transition rate r p in the stationary state assuming that the number of planar 
membranes is much larger than the number of vesicles. This can be done by assuming a two-state 
dynamics in the presence of a barrier. In this case the net current 



j = —a < 27TK, 



sin(0) + - cos , 
y ' 2 \2 



P + ^bT^}, (.11) 
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obtained from the Fokker-Planck equation (1130 . is a constant in #-space, in the stationary case. 
Then, the transition rate r p is defined by r p = j/P p , where P p is the total number of planar 
membranes that can be calculated by integrating the stationary distribution function P s (9) from 



r to 9 = it 



36]. 



The explicit expression for j and P p can be obtained by expanding in a Taylor series up to 
second order in 9 the free energy potential F about its maximum at 9* and its local minimum at 
9 = ir. This procedure yields the approximate expressions 

F max (0) ^ ^ (16 + 5 2 ) - ^(16 - S")(9 - 9*) 2 , (15) 
4 'lb 

F w (0) ~2TrK b 8 + ^(4-5)(9-Tr) 2 (16) 

To calculate j and P p one then uses Eqs. (fT5|) and (fl~6|) . respectively. Now, the transition rate 
r p from planar membranes to spherical vesicles is 



rp = _ S)(4 + 6) 1 / 2 e -™*(*-4) a /« B r (17) 

This expression is valid as long as the energy barrier is larger than the thermal energy. According 



to our expression this imposes the condition 5 < 4 — y 4 ^^ . When this condition is not satisfied, 
the vesicle formation must be analyzed as a transport process. In Figure [3^ we show r p as a 
function of 5 for different values of K^jhsl '• 

In order to calculate the rate r v of the inverse process when the initial condition is such that 
the number of vesicles is much larger than the number of planar membranes, we first approximate 
the free energy around the local minimum at 9 = 0, obtaining 

F o W^^(l6 + 5 2 )-vrH^-0 2 . (18) 

Note that the quadratic term in the approximation of free energy at 9 = is negative. Thus, 
when evaluating the number of vesicles P v around this minimum we obtain 

^ (19) 
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a=2 
a-15 
a-25 
a=50 




Limit of validity of the ' 
■ ation dynamics approximation 



a=2 




FIG. 3: a) Dimcnsionless transition rate r p as a function of S for different values of a = Kb/kgT = 
2, 15, 25, 50 at room temperature. The vertical lines indicate the limit of validity of the activation dynamics 
approximation, b) Dimcnsionless transition rate r v as a function of 6 for different values of a = 2, 5, 15, 25 
at room temperature. 

where a v d9 is the number of membranes between and d6. In this case, the integration of the 
Boltzmann factor, exp(— F$(9) /ksT), must be evaluated between the minimum and the position 
of the maximum at 9*. For simplicity sake, we have approximated 9* up to first order in 5; that 
is, 9* ~ 5/2. After calculating the integral our estimation for the transition rate r v from spherical 
vesicles to planar membranes is 



In Figure [3l we show the behavior of r v as a function of 5 for different values of a = 2iTKb/kBT '. 
For a given value of a, the transition rate grows as the energy barrier decreases (8 increases). For 
a given value of 5, the transition rate depends on the relative value of the bending energy with 
respect to the thermal energy, a. For constant temperature T, decreasing the bending modulus Kb 
favours the formation of vesicles associated with increasing values of r v . It is interesting to notice 
that the Arrhenius law is a consequence of two ingredients: i) the assumption of the linear law 
Eq. (|12p which implies a process not too far from equilibrium, and ii) the fact that the minimum 
of the potential is an extremum, see Eq. (|17p . This second condition is not fulfilled in the case of 
Eq. (I20p which clearly is not an Arrhenius type law. Non- Arrhenius behaviors can also emerge as 




(20) 




FIG. 4: Angle 6 as a function of the dimensionlcss time t = 2irKbat for the following values of 5: 15, 25 and 
60 obtained by solving numerically Eq. (|2ip. The solid circles represent an analytical solution given in the 
text for 5 = 15. The inset shows a log-log plot of vesicle formation time r c versus 5 for 5 = 10, 15, 25, 60, 
400. The slope of the straight line is ~ —1.07. 



a consequence of other mechanisms, see for example 



37|, 



38]. 



B. Vesicle formation as a transport process 

For values of 5 larger than 4, the absolute minimum of the free energy occurs at 8 = (corre- 
sponding to spherical vesicles) without the presence of energy barriers. Therefore, in this case the 
dynamics must be analyzed transport process. 

To do this, one may neglect thermal fluctuations so that the distribution function can be 
approximated by a Dirac delta function P(0,t) = 5[8 — 9{t)], [39]]. In this case, after multiplying 
by 8 and integrating over all 0-space, Eq. (]13p reduces to the dynamical equation 

d 5 (8 

— air) = — cos - 
dr y ' 2 V 2 

where we have defined the dimensionless time r = I'KQ.K^t. Since this equation cannot be solved 
analytically, we have solved it numerically by using a Runge-Kutta method. The solutions (open 
symbols with lines) as a function of r for three different values of 5 are shown in Fig. SI As initial 
condition we used 6(0) = 3.14 since this value represents a nearly planar membrane but with a 
small perturbation that permits the membrane to evolve to its equilibrium state (closed vesicle). 

As it is clear from the figure, during most of the evolution time the value of 8 is close to ix. 
Therefore, one may approximate Eq. (|2ip around 8 = it to first order in its series expansion. This 
can be solved with the same initial condition leading to 8(r) ~ it — 0.0016 exp[(5 — 4)r/4]. This 
solution is represented in Fig. [J] by the solid circles. As can be seen, it constitutes an excellent 
approximation at all times. 



+ sin (8) 



(21) 
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FIG. 5: Probability distribution P as a function of angle and dimcnsionless time r obtained by numerically 
solving Eq. (| 1 3|) with an initial condition given by a Gaussian distribution centered at 9 = it. At short 
times (r < 1) diffusion dominates spreading the distribution. For times r > 1 the drifting force dominates 
and the distribution becomes narrow close to 9 = 0. 

For large values of 5, this approximation allows us to estimate the characteristic vesicle formation 
time t c , given by t c = 2/(irKi,a 5) = 2/(^\/ttA a 7^ . The inset of Fig. H] shows that the dependence 
of log(r c ) as a function of log(<5) is linear with a slope close to —1. 

We now note that the Onsager coefficient a is a mobility that depends on the dynamic viscosity 
of the solvent r] and the membrane area A. Since a -1 has dimensions of energy multiplied by time, 
it must depend on the combination a -1 ~ rjA 3 ^ 2 , yielding 



This relation predicts that the formation time is given by the ratio between the friction force 
~ r]A exerted by the solvent on the membrane by the linear tension force 7 at the contour of the 
membrane. 

For typical values of vesicle area A ~ 10 5 nm 2 , the dynamic viscosity of the solvent 77 ~ 10~ 3 Pas 
and the linear tension 7 ~ 1 — 2/cgT '/nm, one obtains that the characteristic formation time is of 
the order t c ~ 1 ms. 

Finally, we have numerically solved Eq. (|13p in order to study the effects of thermal fluctua- 
tions. We have taken 5 = 5 and k^T ' /2nKb = 1. As an initial condition we have taken a Gaussian 
distribution centered at 6 = tt. At short times (r < 1) diffusion dominates spreading the distri- 
bution whereas for times r > 1 the drifting force dominates and the distribution becomes narrow 
close to 8 = 0. These results are shown in Fig. [5j 



2 riA 




(22) 



IV. CONCLUSIONS 



In this article, we proposed a free energy expression accounting for the formation of spherical 
vesicles from planar membranes. This energy depends on a single state variable and contains two 
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physical parameters related to the membrane rigidity and to the edge tension. The equilibrium 
properties of this energy depend on the ratio, 5 = (jrA) 1 ' 2 ^/^, between the contour energy and 
the Helfrich curvature energy. When 5 < 4 the free energy presents a barrier which disappears for 
larger values of 8. 

Using mesoscopic nonequilibrium thermodynamics rules and the equilibrium information, we 
have derived a Fokker-Planck equation for the probability distribution describing the dynamics of 
vesicle formation. 

Two cases have been analyzed: i) Formation in the presence of barriers (5 < 4) and ii) formation 
as a transport process (<5 > 4). In the first case we have derived expressions for the transition rates 
of formation of vesicles from planar membranes (r p ) and viceversa (r v ). Our expression for r p 
follows an Arrhenius law [see, Eq. (|17p ] and is an increasing function of 5. The rate of the inverse 
process r v has an unusual dependence on temperature [see, Eq. ([20]) ] due to the fact that the free 
energy minimum at 8 = (vesicles) is not an extremum. We have found that r v /r p is orders of 
magnitude smaller than one thus implying that the unwrapping of the spherical vesicles is a very 
improbable process even in the case when the free energy favours it. 

In the second case, the free energy minimum always corresponds to spherical vesicles and can be 
analyzed by using a deterministic equation for the angle as a function of time after neglecting the 
effects of thermal fluctuations. A simple analytical expression that is an excellent approximation 
of the numerical solution allows us to estimate the characteristic vesicle formation time t c which is 
proportional to the membrane area and the viscosity of the solvent, and inversely proportional to 
the edge tension 7, see Eq. (I22p . For typical values of phospholipid membranes of linear dimensions 
~ lOOnm, we have obtained that t c ~ 1ms. 

These results suggest that in typical experiments involving video-based measurements, the dy- 
namics of vesicles formation is dominated by a transport process whereas for numerical simulations, 
in which the studied system is small, the presence of energy barriers could be relevant to the dy- 
namics. 

It could be interesting to test the present model by performing single vesicle experiments and 
simulations in which the detailed evolution in time of the membrane edge can be followed, so that 
the characteristic vesicle formation time can be obtained. 

The proposed model could be useful in the understanding of the mechanisms of phospholipid 
vesicle formation widely used as model experimental systems to study the thermoelastic properties 
of cellular membranes. 
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